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Abstract 

We simulate the self-propulsion of devices in a fluid in the regime of low 
Reynolds numbers. Each device consists of three bodies (spheres or cap- 
sules) connected with two damped harmonic springs. Sinusoidal driving 
forces compress the springs which are resolved within a rigid body physics 
engine. The latter is consistently coupled to a 3D lattice Boltzmann frame- 
work for the fluid dynamics. In simulations of three-sphere devices, we find 
that the propulsion velocity agrees well with theoretical predictions. In sim- 
ulations where some or all spheres are replaced by capsules, we find that the 
asymmetry of the design strongly affects the propelling efficiency. 

Keywords: Stokes flow, self-propelled microorganism, lattice Boltzmann 
method, numerical simulation 



1. Introduction 

Engineered micro-devices, developed in such a way that they are able to 
move alone through a fluid and, simultaneously, emit a signal, can be of cru- 
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cial use in various fields of science. A few engineered self-propulsive devices 
have been developed over the last decade. In one approach, the device con- 
sisted of magnetically actuated superparamagnetic particles linked together 
by DNA double strands [8]. In another approach, miniature semiconductor 
diodes floating in water were used. When powered by an external alternat- 
ing electric field, the voltage induced between the diode electrodes changed, 
which resulted in particle-localized electro-osmotic flow. In the latter case, 
the devices could respond and emit light, thus providing a signal that could 
be used as a carrier [5j . A third example is a swimmer based on the creation 
of a traveling wave along a piezoelectric layered beam divided into several 
segments with a voltage with the same frequency but different phases and 
amplitudes applied to each segment [21j. 

Propulsion through a fluid induces the flow of the surrounding fluid, which 
in turn affects the propulsion of the device. Most generally the fluid motion 
is described by the Navier-Stokes equation 

+ (u • V)u = — Vp + ^ V 2 u + /. (1) 
at p p 

Here u is the velocity, p is the density, p is the pressure, p is the dynamic 
viscosity of the fluid, and / is the applied force in the fluid. The first term on 
the left hand side describes the acceleration of the fluid, whereas the second 
term accounts for non-linear inertial effects that may give rise to effects such 
as turbulence. On the right hand side the first term is the driving pressure 
gradient, while the second term takes into account the viscous dissipation. 

In the case of low Reynolds numbers, the inertia terms on the left side of 
equation ([T]) can be dropped and one is left with 

= --Vp+^V 2 £+/, (2) 
p p 

which together with the incompressibility equation V • u = describes the so- 
called Stokes flow. The main characteristics of the Stokes flow emerge from 
the domination of viscous forces. Consequently, the flow is always laminar 
(no turbulence or vortex shedding) and characterized by a small momentum. 
Since the flow is proportional to the forces applied, linear superposition of so- 
lutions is valid. Furthermore, the Stokes flow is instantaneous, which means 
that the time-development is given solely by the effects of the boundaries. 
Finally, there is no coasting, and the flow is time-reversible, which has im- 
portant implications for swimming strategies at low Reynolds numbers [36J. 
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The propulsion strategy in this regime must involve a non-time-reversible 
motion, thus involving more than one degree of freedom [32]. Such a re- 
quirement is not limiting if a large number of degrees of freedom is available. 
However, if one tries to design a self-propelled device (e.g. a swimmer) with 
only a few degrees of freedom, a balance between simplicity and functionality 
has to be found. 

A number of analytic and numerical studies have been performed to un- 
derstand the behavior of a single or a couple of swimmers under various 
conditions p~3|,[22]. However, most of these types of approaches are limited 
to simple geometries of swimmers and fail when numerous swimmers are in- 
volved leading to a collective behavior. In addition to the effective treatment 
of many body interactions, the regimes in which inertia starts to play a role, 
as well as the problem of motion in confined geometry or transport in tur- 
bulent flow, are inaccessible. In these circumstances extensive simulations 
become essential, but difficult to perform efficiently [3T1 138] . 



Figure 1: A three-sphere swimming device with two translational degrees of freedom. 

In order to meet these challenges, we augmented our already existing mas- 
sively parallel lattice Boltzmann simulation framework waLBerla, widely 
applicable Lattice Boltzmann solver from Erlangen [10J, with swimmers. The 
motion of the latter is included into waLBerla by coupling it with the pe 
rigid multibody physics engine, a framework for the simulation of rigid, fully 
resolved bodies of arbitrary shape [16]. This allows us to resolve both, the 
driven motion of a swimmer within the fluid, and the induced motion of 
the fluid in a consistent manner. For the swimming device we choose the 
simplest possible design, consisting of three rigid bodies connected with two 
springs (Figure [j]). This design has been studied extensively in the three 
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spheres geometry by analytical methods fT2] . 

This paper starts with the elucidation of the background of the numeri- 
cal methods, which form the basis of our simulation, in Section [2j Here, the 
lattice Boltzmann method (LBM) of the waLBerla framework, the me- 
chanics of the pe framework and, finally, the coupling of both with special 
focus on the swimmers are briefly explained. In Section [3| our new im- 
plementation approach for the integration of the swimmers is considered in 
detail, including the cycling strategy. Moreover, our design variants are pre- 
sented. Afterwards, in Section [IJ a quantitative validation of the simulation 
results of the pe engine is conducted. In Section [5j validation criteria for the 
swimmer are defined, and asymmetric and symmetric designs are simulated 
and analyzed in detail to demonstrate the capability of our extension of the 
framework. Finally, in Section [6] a conclusion of the achieved results is given, 
and propositions for future work are made. 

2. Numerical methods 

2.1. The lattice Boltzmann method 

The LBM can be used as an alternative to the classical Navier-Stokes 
(NS) solvers to simulate fluid flows [£3, [35] The fluid is modeled as a set 
of imaginary, discretized particles positioned on an equidistant grid of lat- 
tice cells. Moreover, the particles are only allowed to move in fixed, pre- 
defined directions given by the finite set of velocity vectors, resulting from 
the discretization of the velocity space. Unlike other methods, which rely 
on the automatic generation of a body fitted mesh during the simulation, 
the waLBerla framework uses the LBM on the same (block-)structured 
mesh, even when the objects are moving. This eliminates the need for dy- 
namic re-meshing, and hence offers a significant advantage with respect to 
performance. 

In this study we use the common three-dimensional D3Q19 velocity phase 
discretization model originally developed by Qian, d'Humieres and Lalle- 
mand [34] with TV = 19 particle distribution functions (PDFs) f a :QxT^ 
[0,1), where C M 3 and T C are the spatial and the time domain, 
respectively. Figure [2] illustrates the 19 directions. The corresponding di- 
mensionless discrete velocity set is denoted by {e a \a — 0, . . . , N — 1}. This 
model has been shown to be both stable and efficient [26] . For the work pre- 
sented in this paper, we adopt a lattice Boltzmann collision scheme proposed 
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Figure 2: The D3Q19 velocity phase discretization model. 

by Bhatnagar, Gross and Krook, called lbgk j3l [34] 

f a (xi + e a At : t + At) = f a (x h t) - -[f a (x h t) - f^\x u t)} , (3) 

T 

where Xi is a cell in the discretized simulation domain, t is the current time 
step, t+ At is the next time step, r is the relaxation time in units of time step 
At (which is set to be 1 here) and fa^ represents the equilibrium distribution. 

The time evolution of the distribution function, given by equation ([3]), is 
usually solved in two steps, known as the collision step and the streaming 
step, respectively: 

U(xi,t) = f a (x h t) - -\f a (xi,t) - f^(xi,t)] , (4) 

T 

f a (x t + e a At, t + At) = f a (x h t) , (5) 

where f a denotes the post-collision state of the distribution function. The col- 
lision step is a local single-time relaxation towards equilibrium. The stream- 
ing step advects all PDFs except for f to their neighboring lattice site de- 
pending on the velocity. Thus, for each time step only the information of the 
next neighbors is needed. 

As a first-order no-slip boundary condition a simple bounce-back scheme 
is used, where distribution functions pointing to a neighboring wall are just 
reflected such that both normal and tangential velocities vanish: 

U(xf,t) = /a (£/,*) , (6) 
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with a representing the index of the opposite direction of a, = — e a , and 
Xf explicitly denoting the fluid cell. 

Due to the locality of the cell updates, the LBM can be implemented 
extremely efficiently (see for instance [30l E7j). For the same reason, the 
parallelization of LBM is comparatively straightforward [20] . It is based on 
a subdomain partitioning that is realized in the waLBerla framework by a 
patch data structure as described in Feichtinger et al. [10J. 

2.2. The rigid body physics engine pe 

Simulating the dynamics of rigid bodies involves both the treatment of 
movement by a discretization of Newton's equations of motion, as well as 
handling collisions between rigid objects. The forces involved may be ex- 
ternal, such as gravity, or between bodies, such as springs. Furthermore, a 
number of velocity constraints are available through which hinges, sliders, 
ball joints and the like are implemented [28]. The collisions between rigid 
bodies are either treated directly by the application of the restitution laws 
in the form of a linear complementarity problem (LCP) [2j or by applying 
Hertzian contact mechanics (as for instance in the DEM approach) [7j. 

In our work, we fully resolve the geometry of the rigid bodies. In con- 
trast to mass-point-based approaches, this for instance enables us to easily 
exchange the geometries of the swimmers in our simulations (see Subsec- 
tion 3.2). The rigid body framework used for this is the pe rigid multibody 
physics engine [16J. Due to its highly flexible, modular and massively parallel 
implementation, the framework allows for a direct selection of the time dis- 
cretization scheme and collision treatment. By this it can easily be adjusted 
to various simulation scenarios. For instance, it has already been successfully 
used to simulate large-scale granular flow scenarios [TT] and can be efficiently 



coupled for simulations of particulate flows [15] (see Subsection 2.3). 

The springs of the pe are modeled as damped harmonic oscillators in 
which the bodies are subject to the harmonic and the damping force. The 
first one is 

Spring = ~ k 5n for all springs n, (7) 

where k is the force constant and Ax Sn = (x Bi — x Bj ) — Iq is the deformation of 
the spring for i ^ j. The deformation is defined as the difference between the 
rest length l of the spring and its current length, which equals the current 
distance between the connected bodies Bi and Bj ) whose current positions are 
denoted by x Bi and x Bj , respectively. In order to account for the dissipation 
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in the spring, a damping force, proportional to the velocity Au Sn at which 
the spring is contracting is implemented as in the following 

^cfamp —jAu Sn for all springs n. (8) 

Here 7 is the damping coefficient and Au Sn = u Bi — u Bj for i ^ j. Here, u B% 
and u Bj are the velocities of the bodies Bi and Bj . This damping force hence 
models the internal friction resulting from the deformation of the spring. 
Consequently, each of the two bodies Bi and Bj connected by a spring S n 
feels the force 



F B i { = —k Ax Sn - 7 Au Sn , (9) 
= k Ax Sn + 7 Au Sn , (10) 

respectively. 

A standard way to classify a damped harmonic oscillator is by its damping 
ratio 

lym k 

If D > 1 the system is overdamped and upon an excitation the system returns 
exponentially to equilibrium without oscillating. If D <C 1 the system is 
underdamped. It will oscillate with a re-normalized frequency decreasing 
the amplitude to zero until the equilibrium state is reached. We choose the 
parameters of the springs of the swimmer in such a way that the latter regime 
is achieved. 

Algorithm [T] gives an overview of the necessary steps within a single time 
step of the pe framework. During each time step, collisions between rigid 
bodies have to be detected and resolved appropriately in order to keep rigid 
bodies from interpenetrating. Therefore, initially, all contacts k between 
pairs of colliding rigid bodies are detected and added to the set of constraints 
c, which already includes the potentially existing velocity constraints. In the 
subsequent constraint resolution step, the acting constraining forces F B \ 
including the forces of the harmonic potentials F B ^^ are determined. More- 
over, the detected collisions are resolved depending on the selected collision 
response algorithm. It is important to note that interpenetration cannot be 
completely prevented due to the discrete time stepping, but small penetra- 
tions can be corrected in the collision response phase. In the final step all 
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Algorithm 1 Rigid Body Time Step 
1: / /Collision detection 
2: for each rigid body Bi do 
3: Detect all contacts k 
4: Add them to set of constraints c 
5: end for 
6: 

7: // Constraint resolution 

8: for each constraint c do 

9: Determine acting constraint forces 

10: end for 

11: 

12: / /Time integration 

13: for each rigid body Bi do 

14: Apply forces Pf \ : update position and velocity 
15: end for 



rigid bodies Bi are moved forward in time according to Newton's laws of 
motion, depending on their current velocity and the total forces F t f| acting 
at the given instance of time. 

In the simulations of swimmers, the size of the time step was chosen such 
that rigid bodies do not interpenetrate. Another limitation on the size of the 
time step originates from the force constant of the harmonic oscillators and 
the frequency of the driving forces. In all cases the time step must be small 
enough to enable all bodies to respond sufficiently quickly to changes in the 
acting forces. 

2.3. Coupling lattice Boltzmann to the pe for swimmers 

The coupling between the lattice Boltzmann flow solver and the rigid body 
dynamics simulation has to be two-way: rigid bodies have to be represented 
as (moving) boundaries in the flow simulation, whereas the flow corresponds 
to hydrodynamic forces acting on the rigid bodies (Figure [3]). We use an 
explicit coupling algorithm, shown in Algorithm [2| which we adopt for the 
swimmers. The additional steps are highlighted. Here, a driving force acts 
on the corresponding bodies according to a carefully defined protocol (see 



Subsection 3.3). 



The first step in the coupled algorithm is the mapping of all rigid bodies 
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Figure 3: Illustration of the two-way coupling of waLBerla fluid solver and pe rigid 
body dynamics solver. 
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(a) Initial setup: The ve- 
locities u of the object 
cells x\y are set to the 
velocity u w (xb) of the 
object. In this exam- 
ple the object only has 
a translational velocity 
component. Fluid cells 
are marked with x f . 




(b) Updated setup: Two fluid 
cells have to be transformed 
to object cells, and for two 
object cells the PDFs have to 
be reconstructed. 



Figure 4: 2D mapping example [13] . 
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Algorithm 2 Coupled LBM-PE solver for the swimmers 
1: for each swimmer S do 
2: Map S to lattice grid 
3: end for 

4: 

5: for each lattice cell do 
6: Stream and collide 
7: end for 

8: 

9: for each surface cell do 

10: Add forces F^j dro from fluid to rigid objects Bi 

11: end for 

12: 

13: for each body Bi of the swimmers do 
14: Add external forces 
15: end for 
16: 

17: / /Time step in the rigid body simulation (see Algorithm [j]) 

18: for each rigid body Bi do 

iq. j? B * — fi B * _ i— fi B * _ i— p B i 

iy - r *o* — r hydro ' r dri ' r c 

20: end for 



onto the lattice Boltzmann grid (see Figure 4(a) for an example). Objects 
are thus represented as flag fields for the flow solver. In our implementation, 
each lattice node with a cell center inside an object is treated as a moving 
boundary. For these cells, we apply the boundary condition 

/a (£/,*) = f a (xf,t) + 6w a p w (x f ,t)ea • u w (x f + e a ,t) , (12) 

which is a variation of the standard no-slip boundary conditions for moving 
walls j39]. Here p w {xf ) t) is the fluid density close to the wall, w a is the weight 
in the LBM depending on the stencil and the direction, and the velocity 
u w (xf + e ai t) of each object cell corresponds to the velocity of the object 
at the cell's position at time t. In this way, the fluid at the surface of an 
object is given the same velocity as the object's surface which depends on 
the rotational and the translational velocities of the object. 

The final part of the first step is to account for the flag changes due to the 
movement of the objects. Here, two cases can occur (see Figure [4(b)] ): the 
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fluid cells Xf can turn into object cells, resulting in a conversion of the cell 
to a moving boundary. In the reverse case, i.e. a boundary cell turns into a 
fluid cell, the missing PDFs have to be reconstructed. In our implementation, 
the missing PDFs are set to the equilibrium distributions fa Q \p, u), where 
the macroscopic velocity u is given by the velocity u w (xf) of the object cell, 
and the density p is computed as an average of the surrounding fluid cells. 
Further details of this procedure have been discussed previously by us in 
Iglberger et al. [18]. There are more formal ways for the initialization of the 
missing distribution functions. One of them is described in the paper by Mei 
et al. [25]. 

During the subsequent stream and collide step, the fluid flow acts through 
hydrodynamic forces on the rigid objects: Fluid particles stream from their 
cells to neighboring cells and, in case they enter a cell occupied by a moving 
object, are reversed, causing a momentum exchange between the fluid and 
the particles. The total hydrodynamic force i^ dro on each object Bi resulting 
from this momentum exchange [39] can be easily evaluated due to the kinetic 
origin of the LBM by 

19 ^ 

^hydro = $ZX^ a [ 2 /a +6w a p w (x f ,t)ea-u w (xf + e cn t) , (13) 

where x^ are all obstacle cells of the object that neighbor at least one fluid 
cell. A comparison of different approaches for the momentum exchange is 
given in Lorenz et al. [24] . 

In the third step, the driving force , resulting from the cycling strategy, 
is added to the previously calculated hydrodynamic force on each body Bi. 
Together with the constraint force this results in a total force 

Fit = F*; d ro + *S + F? (14) 

on each object. 

The fourth and final step of the coupling algorithm consists of determining 
the rigid body movements of the objects due to the influence of F^ Q { in the 
rigid body framework. This results in a position change of the objects, which 
then are mapped again to the LBM grid in the next time step. A validation 
of this method is described in Binder et al. [4J, and Iglberger et al. [18]. 
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3. Integration of the swimming device 

The modeling of a three-sphere swimmer always involves a similar physi- 
cal setup. On the one hand, the connections between the objects need to be 
modeled. Here, stiff rods [271 E3] or springs [H] are most commonly used. 
On the other hand, the cycling strategy, responsible for the characteristic 
movement of the swimmer, needs to be defined. 

3.1. Elementary design 




Figure 5: A three-sphere swimming device. The three spheres have the same mass m sp h 
and radius r sp h. Similarly, the two damped harmonic oscillators have identical force con- 
stant fc, damping parameter 7, and rest length Iq. 



Our elementary design (Figure [5]), inspired by the analytical modeling in 
Golestanian and Ajdari [12], and in Felderhof jTTj, consists of three spheres 
(£>i, £>2, £>3) with identical radius r sp h and mass m sp h. The connections be- 
tween the three rigid objects of our swimmer are realized by two damped 
harmonic springs Si and S2 with equal force constant k and damping pa- 
rameter 7. The forces are applied on the spheres along the main axis of the 
swimmer, in our case the z-axis, resulting in a translation of the swimmer in 
this direction. 

In accordance with the equations ^ and (10), the different objects feel 
the force 

F^ ci = -k Ax Sl - 7 Au Sl , (15) 
= k (Ax Sl + Ax^ 2 ) + 7 (A^ 1 + A^ 2 ) , (16) 
= -k Ax S2 - 7 A^ 2 , (17) 

respectively. 
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Since the bodies of the swimmer do not collide with each other, the total 



force acting on each body Bi of the swimmer, given in equation (14), can be 
rewritten as 



r tot — r hydro ' r dri 



(18) 



3.2. Alternative designs 

Since we want to investigate the effect of different geometries on the swim- 
ming behavior of our self-propelled device, we will replace the spheres of the 
elementary three-sphere swimmer design by capsules. For example, Figure [6] 
depicts the characteristic parameters of a three-capsule swimming device. 
The comparably smooth edges of a capsule, and thus its resulting smooth 
behavior in a fluid, substantiate our choice of this particular geometry. Fig- 
ure [7] shows our assembly variants. 




Figure 6: A three-capsule swimming device. The three capsules have the same mass m cap , 
length / cap and radius r cap . Similarly, the two damped harmonic oscillators have identical 
force constant /c, damping parameter 7, and rest length Iq. 



3.3. Cycling strategy 

Purcell's Scallop theorem [32] necessitates the motion of a swimmer to 
be non-time-reversible at low Reynolds numbers. Additionally, the total 
applied force on a swimmer should vanish over one cycle, which means that 
the displacement of the center of mass of a swimmer over one cycle should 
be zero in the absence of a fluid. 

In contrast to the common approach of imposing known velocities on the 
constituent bodies of a swimmer, used, e.g., by Earl et al. [9j and Najafi and 
Golestanian [27], our driving protocol imposes known forces on the objects. 
These forces, only applied along the main axis of the swimmer on the center 
of mass of each body, are given by 
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(19) 

(20) 
(21) 



Here a is the amplitude, uj is the driving frequency, T is the oscillation 
period, and is the phase shift, kept constant at a value of tp = T/4. In 



equation (21), we apply the negative sum of the forces of the two outer bodies 



B 2 and B 3 on the middle body B\ in order to ensure that the net driving 
force acting on the system at each instant of time is zero. Figure [8] illustrates 
these forces. 
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Figure 8: Fragments of the z-components of the driving forces on all three bodies. 

Applying this force protocol to our swimmer, we end up with our cycling 
strategy illustrated in Figure |9| Starting from the rest position, we have to 
delay the force on body £> 3 by a fourth of the pulse length. Step (i) depicts the 
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(a) Theoretical positions prescribed (b) Position plot from the simulation. The z- 
by the cycling strategy. positions are given in lattice cells. 

Figure 9: Cycling strategy of the non-time-reversible motion, exemplified on the three- 
sphere swimmer. For both damped harmonic potentials, l m [ n is the minimum length of 
one arm, l is the rest length, and / max is the maximal extended armlength. A is the 
distance covered by the swimmer after one cycle. 



initial position of our swimming device. The two "arms " , i.e. the connecting 
harmonic oscillators S\ and S 2 , are at their rest length / - In the following, 
li and l 2 denote the current lengths of Si and £2, respectively. From step 
(i) to step (ii), we first apply the driving force on body B 2 and F^ r \ 
on body Bi. As shown in Figure [8j and also apparent from equation ( [l9| ), 
F^? starts in the negative direction until it reaches its negative maximum 
at which point the oscillator Si obtains its maximum length / max . S 2 is also 
influenced by the force and, therefore, also gets extended to the length 
lo < h < ^max- From step (ii) onwards, F^? increases and we start to exert 
the positive driving force F^\ on body £> 3 . We reach the intermediate step 
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(iii), where both oscillators Si and S 2 have the length l < li = l 2 < Z max . In 
the transition from step (iii) to step (iv), F^? goes to 0, i.e. Si has decreasing 
length l < li < / max . Moreover, F^ Y f reaches its positive maximum, and thus 
S2 obtains its maximum length / max . With F^? still increasing and 
starting to decrease, we attain another intermediate state in step (v). Here, 
the length of Si is / m i n < li < Zq, and the length of S2 is lo < h < 'max- 
In step (vi), F^? reaches its positive maximum. As a result, Si takes its 
minimum length / min . Furthermore, F^ Y f goes to 0. Accordingly, Si relaxes 
to the length Z min < li < l . Now, we start to decrease F^? and also decrease 
Ffaf further. Passing by the temporary step (vii), where both oscillators Si 
and S 2 are at the length lo < h = k < Z m in, we pass into step (viii). Here, 
F^ Y f reaches its negative maximum, resulting in £2 obtaining its minimum 
length Z min . Additionally, F^? goes to and therefore Si relaxes to the length 
Z m in < h < lo- While F^? continues to grow in the negative direction and F^f 
decreases, our swimmer moves forward to another transitional state (ix), for 
which the conditions Z < '1 < 'max and / mm < l 2 < lo hold for the oscillators' 
lengths. Finally, F^ 2 reaches its negative maximum again, and F^ Y f becomes 
0, resulting in a length of l < h < 'max for the oscillator 5 2 . On reaching 
the state (x) (which is equal to state (ii)), one swimming cycle is completed. 
Subsequently, our swimmer can begin another cycle of motion. 



4. Validation of the simulation of swimmers with the pe 

In order to quantitatively validate the application of springs in the pe 
physics engine, we perform simulations of the motion of the three-sphere 



swimmer (Figure 7(a) ) in vacuum, and compare these results with the anal- 
ogous analytic model. The basis of that model is the Lagrangian of a non- 
dissipative assembly 



^ _ m xj+xj+xj _ £ (^i-^2-r ) 2 +(x3-xi-r ) 2 j 



(22) 



Here, Xi (i = 1,2,3) refers to the derivative of the positions of the three 
spheres with respect to time and l denotes the rest length of the springs. The 
first term in the Lagrangian consists of the kinetic energy of the three spheres, 
the second gives the energy stored in the two springs due to their deformation, 



and the remaining terms account for the driving forces (equations ( 19)— (21 )). 
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As discussed previously (equation (|8])), the damping forces in the simu- 
lation are proportional to the velocity of the spring deformation, with the 
factor of proportionality being 7. Taking the damping into account, the 
equation of motion for each sphere becomes 

I|=i + ^ fM! ^' 2 ' 3 ' (23) 

This set of three differential equations can be solved analytically assuming 
appropriate initial conditions. 




pe simulation 

Analytical 

calculation 



Time step [10 4 s] t 1 t 2 t 3 



Figure 10: Plot of the z-positions of the three spheres vs. time (weak damping) of one of 
our simulations. The values of the parameters used are: m = 5.44 x 10 -13 kg, k = 1.72965 
kg/s 2 , (J ) z = 25 lattice cells (with discretization Ax = 10" 6 m), 7 = 1.57237 x 10" 7 kg/s, 
uj = 296057.86703 1/s, a = 1.0 x 10~ 5 kgm/s 2 , (x 2 (0)) z = 25 lattice cells, (x 1 (0)) z = 50 
lattice cells, (^3(0))^ = 75 lattice cells, £2(0) = m/s, xi(0) = m/s, 2:3(0) = m/s. 
The simulation is run for 196812 time steps, while the driving force on the left body is 
switched off at t\ = 140580 time steps, and on the right and hence also on the middle body 
at t 2 = 147609 time steps. The springs obtain their starting configuration at t% « 179000 
time steps. 
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Figure [TO] shows both the simulation and the analytically calculated re- 
sults for a typical variation of the sphere positions with time in vacuum, 
in the case of weak damping. The used parameter values are stated in the 
caption. This particular setup of the damped harmonic oscillators and the 
driving forces is equal to the setup for the fluid simulation in Section [5j For 
this number of time steps the lattice Boltzmann algorithm is stable and we 
also end up in the Stokes regime. As expected, each sphere performs an 
oscillatory motion in the steady state with its frequency being the same as 



the driving frequency. Moreover, Figure [10] indicates the force neutrality of 
our total cycle. After we switch off the driving force on the left body at 
t\ = 140580 time steps and on the right body and hence also on the middle 
body at t 2 — 147609 time steps (which equal 3/4 of the total time steps), the 
springs in our swimmer start to revert to their rest lengths, and the swimmer 
soon obtains its starting configuration at t 3 « 179000 time steps, causing the 
flat ends of the position curves. 
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Figure 11: Amplitude of oscillation A scaled by amplitude of the driving forces a (both 
given on the lattice) on the example of sphere 2 vs. uuq/uj, with uo fixed at a value of 
296057.86703 1/s, for very weak damping (D = 0.07). 



Figure [TT] shows the graph of the amplitude of the oscillation A of the 
sphere 2 (scaled by the amplitude of the driving forces a) as a function of 
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the natural spring frequency ujq = ^k/m (scaled by constant a;), with the 
fixed characteristic damping ratio D = 0.07 (equation (11)). 

For this graph, the mass of each sphere was fixed at m — 5.44 x 10 -13 kg, 
while the stiffness k of the spring was varied between the values of 4.76813 x 
10 -3 kg/s 2 and 0.181309 kg/s 2 . Accordingly, the damping 7 lies within the 
bounds of 2.27017 x 10" 9 kg/s and 4.414 x 10" 8 kg/s. Both ranges yield a 
damping ratio of D — 0.07 in each run. In order to guarantee the achievement 
of a steady state in all of the simulations, we perform 10 6 time steps. 



Both Figure [10] and [TT] demonstrate excellent agreement between the sim- 
ulation results and the analytically obtained results. In Figure [TTl the 38 



selected simulation scenarios match the analytically obtained curve with the 
two characteristic maxima at ~ uj and ~ uj/\f?> for not less than 10 
digits. 



5. Simulation of swimmers in a fluid 

In the following we perform simulations of swimmers in a highly viscous 
fluid. Apart from exploring the parameter space in order to achieve simula- 
tions in a physically-stable regime, we also show that the shape of a swimmer 
has a considerable effect on its swimming velocity. 

5.1. Parameters of simulations 

The setup of the domain for the design parameter study is a cuboidal 
channel with (x x y x z) = (100 x 100 x 200) lattice cells, each of which is 
a cube with a side length of 10 -6 m. The z-axis corresponds to the axis of 
movement of the device. The highly viscous fluid has a kinematic viscosity 
of v — 7.36 x 10 -5 m 2 /s and a density of p = 1.36 kg/m 3 , these being 
typical values at room temperature for honey. In order to make the different 
designs comparable, the mass of all spheres and capsules has been set to 
^sph = ^cap = 5.44 x 10 -13 kg. The spheres in the designs (a)-(e), (g)-(j) 
and (l)-(o) all have the radius r sp h = 4 x 10 -6 m. The capsules used in 
the designs (b)-(p) all have the radius r cap = 4 x 10 -6 m and the length 
^cap = 8 x 10 -6 m. In order to observe the influence of the distances to the 
walls of the simulation box, we also consider a three-sphere swimmer (q) with 
larger radii r sp hbi g = 8 x 10 -6 m. 

At t — 0, the swimmer is centered in the channel, with the center of 
the middle body being at the center of the channel with respect to all three 
dimensions. As the time-steps increase, the z-positions of the different bodies 
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of the swimmer change, while the x- and the y-positions do not. Also, at 
£ = 0, the springs are at their rest length of 8 lattice cells (= 8 x 10 -6 m) in 
the designs (a)-(f) and (l)-(q). In order to make the respective armlengths 
of the corresponding designs (b)-(f) and (g)-(k) equal, the springs in the 
designs (g)-(k) have a rest length of 16 lattice cells (= 16 x 10 -6 m) between 
perpendicular capsules, and those between a sphere and a perpendicular 
capsule have a rest length of 12 lattice cells (= 12 x 10 -6 m). The other 
springs in those designs have a rest length of 8 lattice cells (= 8 x 10 -6 m). All 
springs are characterized by a stiffness of k = 1.72965 kg/s 2 and a damping 
constant of 7 = 1.57237 x 10 -7 kg/s. This results in a damping ratio of 
D = 0.07 <C 1, i.e. our system is in the weak damping regime. 

The simulations run for 196812 time steps, performing five swimming 
cycles with an oscillation period of T = 28116 time steps. The sinusoidal 
driving forces on the left and the middle body are applied at T = 0, whereas 
the force on the right body is implemented with a delay of = T/4 = 7029 
time steps. The amplitude of the driving force on the right and the left body 
is a = 1.0 x 10 -5 kgm/s 2 . The force on the middle body is at each instant 
the negative of the sum of the forces on the two other bodies, so that the 
swimmer has no net external driving force at any time step. This means that 
the amplitude of the driving force on the middle body is a for £ < 7029 and 
£ > 140580 time steps, and y/2a otherwise. All simulations are terminated 
by switching off the driving at t\ = 140580 time steps for the left body and 
£2 = 147609 time steps (which equals 3/4 of the total simulation time) for all 
other bodies and the system is allowed to relax on its own for the remaining 
time. It should be noted that the same driving forces are applied on the 
bodies in each simulation, and any differences in their motion are due to 
their different shapes. 

5.2. Stability 

We show that our choice of parameters results in our swimmers being in 
the regime of low Reynolds numbers [32j . The latter is given by 

Re = ^ < 1, (24) 

with the kinematic viscosity of the fluid v ) the characteristic velocity u and 
the characteristic length-scale /. 

In order to simulate our swimmer in a stable regime, we define two dif- 
ferent Reynolds numbers. The first is the Reynolds number of the whole 
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Design 


Rest 
Length 


Total 
Width 


#e swim [10" 3 ] 


Re B * 


Re B - 


Re B * 


(a) 




40 


8 


6.51 


0.066 


0.059 


0.078 


(b) 


48 


8 


6.32 


0.128 


0.047 


0.130 


(c) 


48 


8 


5.75 


0.053 


0.123 


0.081 


(d) 


48 


8 


6.14 


0.074 


0.060 


0.121 


(e) m ^ 


56 


8 


5.72 


0.143 


0.047 


0.130 




64 


8 


5.51 


0.119 


0.096 


0.136 


(g) 




44 


16 


5.07 


0.064 


0.041 


0.086 


(h) 




48 


16 


4.93 


0.047 


0.065 


0.083 


(i) 


44 


16 


4.93 


0.077 


0.060 


0.053 


(j) 


48 


16 


3.64 


0.074 


0.042 


0.059 




56 


16 


3.72 


0.056 


0.044 


0.063 


(1) 




40 


16 


4.86 


0.062 


0.041 


0.085 


(m) 




40 


16 


4.54 


0.045 


0.063 


0.078 


(n) 




40 


16 


4.64 


0.078 


0.058 


0.050 


(o) 




40 


16 


3.64 


0.074 


0.041 


0.054 


(P) 




40 


16 


3.53 


0.051 


0.041 


0.053 


(q) e 


64 


16 


2.33 


0.080 


0.061 


0.083 



Table 1: Parameters of the different swimmers. Re 1 denotes the Reynolds number of 
the ith body in each swimmer, while Re sw i m indicates the Reynolds number of the whole 
swimmer. 
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swimming device (Table [1} fourth column), defined as 

i?e swim = <C 1 (25) 

with respect to the characteristic velocity of the swimmer 2x SW im- The latter is 
calculated from the movement A of the central body in the swimmer during 
one swimming cycle (here chosen to be the last full cycle), divided by the 
time-period T, i.e. ^ sw i m = A/T. The characteristic length-scale / swmi of the 
different swimming devices is the rest length in the direction of motion for 
all the different swimmer designs (Table [TJ second column) . 

Additionally, we verify the Reynolds number not only for the whole swim- 
mer but also for each of the bodies (Table [TJ third, fourth and fifth column), 
which may move with considerably higher velocities than the swimming de- 
vice. Here, the Reynolds number of the body Bi is 

n,Bi ]Bi 

Re Bi = « 1, (26) 

and u Bi is the maximum velocity of the body Bi and l Bi is its length in the 
direction of motion. As Table [T] confirms, all the Reynolds numbers satisfy 
the i?e«l condition. 

5.3. Swimming velocities 

In the following, we investigate the performance of each swimmer as a 
function of its design geometry, as also detailed in Pickl et al. [29]. Specif- 
ically, we study the dependence of the amplitudes of the oscillations of the 
bodies in each swimmer on their shapes, and relate this to the overall swim- 
ming velocity. 

5.3.1. The three- sphere swimmers 

First we look at the two swimmers which consist solely of spheres. Table [2] 
shows the velocity of each swimmer obtained from the simulations (u SW i m ). 
We find that the swimmer consisting of smaller spheres moves significantly 
faster than the one with larger spheres. Figure [12] shows the plot of the 
simulation trajectories of each body in the two swimmers. 

In previous work, this design has been analytically modeled by Golesta- 
nian and Ajdari [12], who predict the velocity of a swimmer (here called 
uqx-) Table [5J third column) if the amplitude of oscillation of its two arms 
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Design 


^swim[lO 


«ga[10" 6 ] 


Error 


Efficiency [1CT 4 ] 


(q) 


2.02 


2.45 


< 18% 


0.88 


(a) 


9.05 


8.98 


< 1% 


10.44 



Table 2: Comparison of swimming velocities from the simulations (u sw i m ) and from the 
prediction of Golestanian and Ajdari («ga), in the case of spheres. All the quantities are 
given in terms of their values on the lattice. 



is known. Here, an arm is defined as the distance between the centers of 
two bodies connected by a spring. In their model, Golestanian and Ajdari 
assume that the deformation of each arm is given by 



hit) = li + di cos(ut + ipi), i = 1, 2 



(27) 



Here ^ is the average length and li(t) the instantaneous length of arm z, 
and di, u and (ft are the amplitude, the frequency and the phase shift of 
the oscillation of arm i. Under these conditions, the formula for the average 
velocity of a three-sphere swimmer is 



^ga = K di d 2 to sin((ft - cp 2 ) 



where K is a geometrical factor given by 



K 



3rir 2 r 3 



2(ri + r 2 + r 3 ) 2 



1 



I 2 



(h + k) 2 . 



(28) 



(29) 



Here, ri, r 2 and r 3 are the radii of the three spheres. It is also assumed that 
Ti <C lj and di <C Tj, for all i and j. 

Since this model assumes a velocity protocol in contrast to our force 
protocol (i.e., we start off with known forces which are applied on the bodies 
whereas the model of Golestanian and Ajdari assumes known deformations 
of the arms, regardless of the forces that cause them), we first extract {di} 
and {cpi} from the simulation trajectories (Table [3]), and use them with the 



known values of {^}, {r^} and uj to find uqa from equation (28). To extract 
{di} and {(ft}, we fit the armlength trajectories for the last complete period 
of motion (from t = 4T to t 

We define the error 
expected uqa value, as 



m u, 



5T) with equation (27). 
swim, when compared with the theoretically- 



Error 



^GA 



x 100. 



(30) 



24 



130 i_i i i i i i i i i | i i i i i i i i i | i i i i i i i i i | i i i i i i i i i | i i i i i i i i i | i i i i 1 1 1 1 1 1 i i i i i i 1 1 1 | i i i i i i i i i_ 




Figure 12: Trajectories of each sphere in the two three-sphere swimmers. The dashed 
curves and the solid curves show the trajectories of the spheres in the large and the small 
swimmer, respectively. The simulation is run for 196812 time steps, while the driving force 
on the left body is switched off at t\ — 140580 time steps, and on the right and hence 
also on the middle body at t<i = 147609 time steps. The springs obtain their starting 
configuration at ts ~ 179000 time steps. 

For the swimmer with large spheres, this error comes out to be 18%. This is 
possibly explained by the effect of the boundaries of the simulation box on the 
spheres of the larger swimmer. The theoretical result assumes infinite fluid, 
whereas the simulations are conducted in a finite box of the dimensions stated 
above. Another possible reason for the disagreement in the two velocities for 
the large-spheres swimmer is that the condition <C lj is not satisfied very 
well. For the swimmer with small spheres, the effect of the finiteness of the 
swimming domain is smaller, and the <C lj condition is also satisfied to 
a greater extent. Consequently, we find that in this case there is excellent 
agreement between uqa and ^ sw i m , to an error of less than 1% (Table [2]). 
We also observe that for the swimmer with small spheres, the amplitude of 
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Design 


K 




d 2 


sin(v?i - Lp 2 ) 


(q) e 


0.0040 


1.59 


2.07 


0.83 


(a) 


0.0045 


2.53 


3.80 


0.94 



Table 3: Parameters for the three-sphere swimmers, uj equals 0.00022. All the quantities 
are given in terms of their values on the lattice. 



oscillation of each arm is larger than that of the corresponding arm for the 
larger swimmer (Table |3j third column)), and that for each swimmer, the 
amplitude of oscillation of the leading arm (rf 2 ) is greater than that of the 
trailing arm (di). 

Table [2] also shows the swimming efficiency of the two swimmers, as de- 
fined by Lighthill [23], 



67T7yr eff ul 
^I^^F^-u^dt 



Efficiency = - r5T ^ - - . (31) 



Here r e ff is an effective hydrodynamic radius of the swimmer, and is approx- 
imated by the sum of the radii of the three spheres pQ . u Bi (t) is found from 
the simulation curves, using the last full swimming cycle (i.e., from t = 4T 
to t = 5T). The integration is also done over this period, in order to work 
with the steady-state values. Table [2] shows that the efficiency of the swim- 
mer with small spheres is greater, as expected due to its significantly greater 
swimming velocity. 

5.3.2. Swimmers with capsules 

The rest of the swimmers contain capsules, and can be separated into 
three distinct families. The first family consists of the swimmers labeled (b)- 
(f) in Table [TJ These contain capsules which are parallel to the swimming 
direction (henceforth referred to as parallel capsules), and have springs of 
a rest length of 8 lattice cells. The second family consists of the swimmers 
labeled (g)-(k) in Table [TJ These contain capsules which are perpendicular to 
the swimming direction (henceforth referred to as perpendicular capsules), 
and have springs of a rest length of 8, 12 or 16 lattice cells. The third 
family consists of the swimmers labeled (l)-(p) in Table [TJ These also contain 
perpendicular capsules, and their springs all have a rest length of 8, making 
these swimmers shorter than the swimmers in the second family. 
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Design 


r-i r\— 6] 
^swim L J- ^ J 


m G a[1U J 


Error 


T7 f r- ; r-i n— 41 

rL/inciency [1U J 




4.78 


4.50 


< 7% 


3.69 


(e) 


5.68 


5.48 


<4% 


5.31 


(c) 


6.66 


6.31 


<6% 


5.86 


(d) 


7.11 


7.16 


< 1% 


7.37 


(b) 


7.32 


7.21 


< 2% 


7.77 



Table 4: Comparison of swimming velocities from the simulations (u sw im) and from the 
prediction of Golestanian and Ajdari (kqa), in the case of the swimmers with parallel 
capsules. All the quantities are given in terms of their values on the lattice. 



Tables [ij [5] and [6] (second column) show the measured ^ sw i m for the swim- 
mers with parallel capsules, with perpendicular capsules and longer springs, 
and with perpendicular capsules and shorter springs, respectively. In all the 
three tables, the swimmers are arranged in the increasing order of t^wim, 
from top-to-bottom. It can be observed that corresponding designs in each 
family occupy the same position in their respective tables. Here, two designs 
are said to be corresponding if one needs to replace the same spheres in a 
three-sphere swimmer to get both the designs, albeit by different kinds of 
capsules (and, possibly, with springs of different lengths). One also observes 
that the speed of swimming decreases with an increase in the number of cap- 
sules. This can be attributed to the fact that the small spheres have greater 
amplitudes of oscillation than the capsules, which means that the lengths of 
arms which contain spheres have greater amplitudes of oscillation than those 
that do not. A greater amplitude of oscillation of the arms results in a higher 
swimming velocity. 

The simulation trajectories of all swimmers with capsules show that ir- 
respective of the number of capsules in the swimmer, each body exhibits a 
smooth sinusoidal movement in the steady state, and one may attempt to 
use the formula of Golestanian and Ajdari (equation ([28])) for the analysis 
of the swimming motion. While this formula was developed assuming the 
bodies to be spherical, one may hypothesize that small deviations from the 
spherical shapes of the bodies will mostly affect the geometrical pre-factor K 
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Design 


^swim[lO 


«ga[10- 6 ] 


Error 


Efficiency [10~ 4 ] 


( k ) ^-%~* 


3.69 


3.95 


< 7% 


2.12 




4.21 


4.66 


< 10% 


3.08 


(h) +~%r* 


5.71 


6.50 


< 13% 


4.24 


(i) 


6.23 


6.75 


<8% 


5.74 


(g) 


6.40 


6.82 


< 7% 


5.97 



Table 5: Comparison of swimming velocities from the simulations (i£ sw im) and from the 
prediction of Golestanian and Ajdari (wqa), m the case of the second family of the swim- 
mers with capsules. All the quantities are given in terms of their values on the lattice. 



Design 




«ga[10- 6 ] 


Error 


Efficiency [10" 4 ] 


(p) 




4.90 


6.70 


< 27% 


3.63 


(o) 


5.05 


6.63 


< 24% 


4.08 


(m) 




6.30 


8.96 


< 30% 


4.76 


(n) 




6.44 


7.74 


< 17% 


5.93 


(1) 




6.75 


7.66 


< 12% 


6.47 



Table 6: Comparison of swimming velocities from the simulations (i£ sw im) and from the 
prediction of Golestanian and Ajdari (i^ga), in the case of the third family of the swimmers 
with capsules. All the quantities are given in terms of their values on the lattice. 
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in equation (28) and not the general expression. We test this hypothesis by 
extracting the amplitudes of the body oscillations and their phase shifts by 
fitting the trajectories with equation (27). The geometric factor is obtained 
by approximating each capsule by a sphere of diameter 12, as that is the 
arithmetic mean of the dimensions of each capsule in the direction of move- 
ment and in a direction perpendicular to it. The mean armlengths are the 
same as in the original swimmers. For ease of reference, the values of all of 
these parameters are given in the appendix (tables A. 7, A. 8 and A.9). These 
tables show that as in the case of the three-sphere swimmers, the amplitude 
of oscillation of the leading arm (d 2 in the tables) is larger than that of the 
trailing arm (d\ in the tables), for each swimmer. 

The second and the third columns of the Tables [4j [5] and [6] show the 
estimated uqa and the error in the measured u sw { m relative to uqa for each 
of the three families. It is interesting to note that the i^ga approximation is 
still in good agreement with the results of the simulations. Especially in the 
first two families the error is about 10% or less. For the third family, i.e. for 
the swimmers with perpendicular capsules and shorter springs, the errors 
increase, to between 10% and 30%, probably because our approximation 
causes the spring lengths to become smaller than or equal to the sphere 
radii. Therefore the conditions necessary for the validity of Golestanian and 
Ajdari's formula are compromised. This is similar to the case of the three- 
sphere swimmer with large spheres, where the calculated and the observed 
velocities of the swimmer also disagreed due to the <C lj condition not 
being satisfied well. 

The tables also show the efficiencies, as defined in equation (31), of the 
swimmers in each family. Here, the hydrodynamic radius A has been cal- 
culated by assuming, as before, that each capsule is replaced by a sphere of 
diameter 12. The results show that in each family, the order of the efficiencies 
accords with the order of u sw i m . This is as expected, since the driving forces 
for each swimmer in all our simulations are exactly the same, and so a more 
efficient swimmer would be expected to move faster than a less efficient one. 

While this comparison has provided some useful insights, in order to 
fully understand our simulation results, the theory of Golestanian and Ajdari 
should be expanded to account for capsules, a task that we hope to undertake 
in the future. 
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6. Conclusions and Future Work 



In this work we have demonstrated the successful integration of self- 
propelled micro-devices into the coupled ^-waLBerla framework. Fur- 
thermore, we have shown the validity of this approach by comparing the 
results of simulations to analytical models. We have taken advantage of the 
flexibility of our framework to simulate symmetric and asymmetric swimmers 
combining capsule and sphere geometries, which has not been possible so far. 

The establishment of this framework is especially important for the study 
of systems that are inaccessible via experimental and analytical methods. In 
the future, we aim to harness its capabilities to address issues such as the 
problem of the three-body swimmer in regimes beyond those dictated by the 
approximations of Felderhof [llj and Golestanian and Ajdari [12]. Another 
interesting application would be to study the behavior of the swimmer in a 
narrow channel and explore the role of the boundary conditions at the wall. 

It would also be interesting to investigate the hydrodynamic interactions 



amongst many swimmers swimming together p~3j[l9]. In Figure 13 we show 
the flow fields averaged over five cycles for the case of a single swimmer and 
also for three swimmers swimming together. The simulation of the three 
swimmers took 60.28 hours of total CPU-time on a single Intel Core i5-680 
3.60GHz. To exclude the effects of the walls of the channel, one needs to 
set up even larger flow fields around the swimmers in this case than in the 
case of a single swimmer. Moreover, our work indicates that a steady state 
is reached after a significantly longer time in the case of many swimmers 
swimming together compared to the case of a single swimmer. 

These factors suggest that in the current setting, the accurate simulation 
of many swimmers would take very long. To reduce this simulation time 
significantly, the large domains of simulation should be processed in a more 
efficient way by parallelizing the framework. Both waLBerla and pe already 
feature massively parallel simulation algorithms for large scale particulate 
flow scenarios [15], but the parallelization of the force generators (i.e. springs) 
and velocity constraints is yet to be implemented. Once this is achieved, 
massive simulations with a large number of swimmers will be accessible, 
allowing us, finally, to address the problem of swarming. 
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(a) 




(b) 



Figure 13: Flow field of one and three three-sphere swimming devices, averaged over five 
swimming cycles. We used the same parameters for the simulation setup as in Section |57[| 
except for the channel, which now is (x x y x z) = (164 x 164 x 240) lattice cells. 
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Design 


K 


di 


d 2 


sm(v?i - (p 2 ) 




0.0030 


2.24 


3.39 


0.89 


(e) ^-«H» 


0.0037 


2.21 


3.63 


0.84 


(c) 


0.0032 


2.59 


3.59 


0.96 


(d) 


0.0041 


2.47 


3.65 


0.87 


(b) 


0.0041 


2.26 


3.79 


0.92 



Table A. 7: The different parameters for the swimmers with parallel capsules, uj equals 
0.00022. All the quantities are given in terms of their values on the lattice. 



Appendix A. Values of different parameters of the swimmers with 
capsules 

Here we present the values of the different parameters of the swimmers 



with capsules, found from the simulation curves and used in equation (28). 
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Design 


K 




d 2 


sin(v?! - <p 2 ) 




0.0030 


2.14 


3.26 


0.85 


(j) 


0.0037 


2.11 


3.62 


0.75 


(h) 


0.0032 


2.66 


3.57 


0.97 


(i) 0-*--% 


0.0041 


2.43 


3.65 


0.84 


(g) 


0.0041 


2.17 


3.80 


0.91 



Table A. 8: The different parameters for the swimmers with perpendicular capsules and 
long springs, uo equals 0.00022. All the quantities are given in terms of their values on the 
lattice. 



Design 


K 


di 


d 2 


sin(<^i - Lp 2 ) 


(p) 




0.0068 


1.92 


2.59 


0.90 


(o) $ • 


0.0058 


2.07 


3.33 


0.76 


(m) 




0.0050 


2.56 


3.24 


0.98 




0.0050 


2.54 


3.39 


0.81 


(1) 




0.0050 


2.00 


3.75 


0.92 



Table A. 9: The different parameters for the swimmers with perpendicular capsules and 
short springs, uo equals 0.00022. All the quantities are given in terms of their values on 
the lattice. 



34 



References 



[1] F. Alouges, A. DeSimone, A. Lefebvre, Optimal strokes for axisymmetric 
micros wimmers, The European Physical Journal E: Soft Matter and 
Biological Physics 28 (2009) 279-284. 

[2] M. Anitescu, F.A. Potra, Formulating dynamic multi-rigid-body con- 
tact problems with friction as solvable linear complementarity problems, 
Nonlinear Dynamics 14 (1997) 231-247. 

[3] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes 
in gases. I. Small amplitude processes in charged and neutral one- 
component systems, Physical Review 94 (1954) 511-525. 

[4] C. Binder, C. Feichtinger, H. Schmid, N. Thiirey, W. Peukert, U. Rude, 
Simulation of the hydrodynamic drag of aggregated particles, Journal 
of Colloid and Interface Science 301 (2006) 155-167. 

[5] S.T. Chang, V.N. Paunov, D.N. Petsev, O.D. Velev, Remotely powered 
self-propelling particles and micropumps based on miniature diodes, Na- 
ture Materials 6 (2007) 235-240. 

[6] S. Chen, CD. Doolen, Lattice Boltzmann method for fluid flows, Annual 
Review of Fluid Mechanics 30 (1998) 329-364. 

[7] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular 
assemblies, Geotechnique 29 (1979) 47-65. 

[8] R. Dreyfus, J. Baudry, M. Roper, M. Fermigier, H.A. Stone, J. Bibette, 
Microscopic artificial swimmers, Nature 437 (2005) 862-865. 

[9] D.J. Earl, CM. Pooley, J.F. Ryder, I. Bredberg, J.M. Yeomans, Model- 
ing microscopic swimmers at low Reynolds number, Journal of Chemical 
Physics 126 (2007) 064703. 

[10] C. Feichtinger, S. Donath, H. Kostler, J. Gotz, U. Rude, WaLBerla: 
HPC software design for computational engineering simulations, Journal 
of Computational Science In Press, Corrected Proof (2011). 

[11] B.U. Felderhof, The swimming of animalcules, Physics of Fluids 18 
(2006) 063101. 



35 



[12] R. Golestanian, A. Ajdari, Analytic results for the three-sphere swimmer 
at low Reynolds number, Physical Review E 77 (2008) 036308. 

[13] R. Golestanian, J. Yeomans, N. Uchida, Hydrodynamic synchronization 
at low Reynolds number, Soft Matter 7 (2011) 3074-3082. 

[14] J. Gotz, K. Iglberger, C. Feichtinger, S. Donath, U. Rude, Coupling 
multibody dynamics and computational fluid dynamics on 8192 proces- 
sor cores, Parallel Computing 36 (2010) 142-141. 

[15] J. Gotz, K. Iglberger, M. Sturmer, U. Rude, Direct numerical simulation 
of particulate flows on 294912 processor cores, in: 2010 ACM/IEEE 
International Conference for High Performance Computing, Networking, 
Storage and Analysis, IEEE, pp. 1-11. 

[16] K. Iglberger, Software Design of a Massively Parallel Rigid Body Frame- 
work, Ph.D. thesis, 2010. 

[17] K. Iglberger, U. Rude, Massively parallel granular flow simulations with 
non-spherical particles, Computer Science - Research and Development 
25 (2010) 105-113. 

[18] K. Iglberger, N. Thiirey, U. Rude, Simulation of moving particles in 3D 
with the lattice Boltzmann method, Computers & Mathematics with 
Applications 55 (2008) 1461-1468. 

[19] A. Kanevsky, M.J. Shelley, A.K. Tornberg, Modeling simple locomotors 
in Stokes flow, Journal of Computational Physics 229 (2010) 958-977. 

[20] C. Korner, T. Pohl, U. Rude, N. Thiirey, T. Zeiser, Parallel lattice 
Boltzmann methods for CFD applications, in: A. Bruaset, A. Tveito 
(Eds.), Numerical Solution of Partial Differential Equations on Parallel 
Computers, volume 51 of Lecture Notes for Computational Science and 
Engineering, Springer Verlag, 2005, pp. 439-465. 

[21] G. Kosa, M. Shoham, M. Zaaroor, Propulsion method for swimming 
microrobots, Robotics, IEEE Transactions on 23 (2007) 137 -150. 

[22] E. Lauga, T.R. Powers, The hydrodynamics of swimming microorgan- 
isms, Reports on Progress in Physics 72 (2009) 096601. 



36 



[23] M.J. Lighthill, On the squirming motion of nearly spherical deformable 
bodies through liquids at very small Reynolds numbers, Communica- 
tions on Pure and Applied Mathematics 5 (1952) 109-118. 

[24] E. Lorenz, A. Caiazzo, A.G. Hoekstra, Corrected momentum exchange 
method for lattice Boltzmann simulations of suspension flow, Physical 
Review E 79 (2009) 036705. 

[25] R. Mei, L.S. Luo, P. Lallemand, D. d'Humieres, Consistent initial condi- 
tions for lattice Boltzmann simulations, Computers & Fluids 35 (2006) 
855 - 862. Proceedings of the First International Conference for Meso- 
scopic Methods in Engineering and Science. 

[26] R. Mei, W. Shyy, D. Yu, L.S. Luo, Lattice Boltzmann method for 3- 
D flows with curved boundary, Journal of Computational Physics 161 
(2002) 680-699. 

[27] A. Najafi, R. Golestanian, Simple swimmer at low Reynolds number: 
Three linked spheres, Physical Review E 69 (2004) 062901. 

[28] K. Pickl, Rigid Body Dynamics: Links and Joints, Master's thesis, 
Computer Science Department 10 (System Simulation), University of 
Erlangen-Niirnberg, 2009. 

[29] K. Pickl, J. Pande, J. Gotz, K. Iglberger, U. Rude, K. Mecke, A. Smith, 
The effect of shape on the efficiency of self-propelled swimmers, submit- 
ted to Physical Review Letters (2011). 

[30] T. Pohl, M. Kowarschik, J. Wilke, K. Iglberger, U. Rude, Optimization 
and profiling of the cache performance of parallel lattice Boltzmann 
codes, Parallel Processing Letters 13 (2003) 549-560. 

[31] CM. Pooley, J.M. Yeomans, Lattice Boltzmann simulation techniques 
for simulating microscopic swimmers, Computer Physics Communica- 
tions 179 (2008) 159 - 164. 

[32] E.M. Purcell, Life at low Reynolds number, American Journal of Physics 
45 (1977) 3-11. 

[33] V. Putz, J.M. Yeomans, Hydrodynamic synchronisation of model mi- 
croswimmers, Journal of Statistical Physics 137 (2009) 1001-1013. 



37 



[34] Y.H. Qian, D. d'Humieres, P. Lallemand, Lattice BGK models for 
Navier-Stokes equation, Europhysics Letters (EPL) 17 (1992) 479-484. 



[35] S. Succi, The Lattice Boltzmann Equation - For Fluid Dynamics and 
Beyond, Clarendon Press, 2001. 

[36] G. Taylor, Analysis of the swimming of microscopic organisms, Pro- 
ceedings of the Royal Society of London. Series A. Mathematical and 
Physical Sciences 209 (1951) 447-461. 

[37] G. Wellein, T. Zeiser, S. Donath, G. Hager, On the single processor 
performance of simple lattice Boltzmann kernels, Computer & Fluids 
35 (2005) 910-919. 

[38] Y.Z. Yang, V. Marceau, G. Gompper, Swarm behavior of self-propelled 
rods and swimming flagella, Physical Review E 82 (2010) 031904. 

[39] D. Yu, R. Mei, L.S. Luo, W. Shyy, Viscous flow computations with the 
method of lattice Boltzmann equation, Progress in Aerospace Sciences 
39 (2003) 329-367. 



38 



